package flexibilityDesign;

import java.math.BigDecimal;

import util.MathUtil;
import util.NullOutputStream;
import ilog.concert.IloLinearNumExpr;
import ilog.concert.IloNumVar;
import ilog.concert.IloRange;
import ilog.cplex.IloCplex;

public class StructureDesign {
	private double lcpm;// c_ij/pi_ij
	private int n;// demand number
	private int m;// resource number
	private int M;// scenario counts

	private double[][] demands;
	private double[] capacities;
	private double[][] c;// link cost
	private double[][] pi;// profit margin

	public StructureDesign() {
	}

	public StructureDesign(int n, int m, int M, double lcpm) {
		this.n = n;
		this.m = m;
		this.M = M;
		this.lcpm = lcpm;
	}

	private void initData() {
		// this.capacities = Support.readCapacities(n, m);
		double cap = BigDecimal.valueOf(n * 100)
				.divide(BigDecimal.valueOf(m), 5, BigDecimal.ROUND_HALF_EVEN)
				.doubleValue();
		this.capacities = new double[m];
		for (int j = 0; j < m; j++)
			this.capacities[j] = cap;
		this.demands = Support.readDemands(M, n);
		this.c = Support.readLinkCost(n, m);
		this.pi = Support.readProfitMargin(n, m);
	}

	public void Lshaped() {
		this.initData();

		double LB = Double.NEGATIVE_INFINITY;
		double UB = Double.MAX_VALUE;
		try {
			long start = System.currentTimeMillis();
			IloCplex master = new IloCplex();
			master.setOut(new NullOutputStream());
			IloNumVar[][] Y = new IloNumVar[n][m];
			for (int i = 0; i < n; i++)
				for (int j = 0; j < m; j++)
					Y[i][j] = master.boolVar();
			IloNumVar theta = master.numVar(-1000000, 1000000);
			int[][] YValue = new int[n][m];// Y value
			double thetaValue = Double.NEGATIVE_INFINITY;
			double masterObjValue = 0;

			IloLinearNumExpr obj = master.linearNumExpr();
			IloLinearNumExpr sumY = master.linearNumExpr();
			for (int i = 0; i < n; i++) {
				for (int j = 0; j < m; j++) {
					obj.addTerm(c[i][j], Y[i][j]);
					sumY.addTerm(1, Y[i][j]);
				}
			}
			obj.addTerm(1.0, theta);
			master.addMinimize(obj);

			master.addLe(sumY, 2 * n);

			int iter = 1;
			int delta = n - n / 2;
			double gap = Double.MAX_VALUE;
			while (true) {
				if (iter > 1) {
					// master.exportModel("D:\\workspace\\CapacityInvestment\\data\\flexibilityDesign\\master"
					// + iter + ".lp");
					IloRange trustRegion = null;
					if (iter < 50) {
						IloLinearNumExpr trustR = master.linearNumExpr();
						int oneC = 0;
						for (int i = 0; i < n; i++) {
							for (int j = 0; j < m; j++) {
								if (YValue[i][j] == 1) {
									trustR.addTerm(-1.0, Y[i][j]);
									oneC++;
								} else
									trustR.addTerm(1.0, Y[i][j]);
							}
						}
						trustRegion = master.addLe(trustR, delta - oneC);
					}

					if (gap > 0.05)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.05);
					else if (gap > 0.03)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.03);
					else if (gap > 0.02)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.02);
					else if (gap > 0.009)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.009);
					else if (gap > 0.007)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.007);
					else if (gap > 0.005)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.005);
					else if (gap > 0.003)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.003);
					else if (gap > 0.001)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.002);
					else
						master.setParam(IloCplex.DoubleParam.EpGap, 0.001);

					double[] initV = new double[m];
					for (int i = 0; i < m; i++)
						initV[i] = 1;
					for (int i = 0; i < n; i++)
						master.addMIPStart(Y[i], initV);

					if (master.solve()) {
						for (int i = 0; i < n; i++)
							for (int j = 0; j < m; j++)
								YValue[i][j] = (int) master.getValue(Y[i][j]);
						masterObjValue = master.getObjValue();
						if (masterObjValue > LB)
							LB = masterObjValue;

						thetaValue = master.getValue(theta);
						// System.out.println("@theta = " + thetaValue);

						if (iter < 50)
							master.getModel().remove(trustRegion);
					} else {
						System.out.println("master status: "
								+ master.getStatus());
					}
				} else {
					// TODO generate initial solution
					for (int i = 0; i < n; i++)
						for (int j = 0; j < m; j++)
							if (j == i)
								YValue[i][j] = 1;
							else
								YValue[i][j] = 0;
				}

				BigDecimal subProblemObjValue = BigDecimal.ZERO;
				BigDecimal et = BigDecimal.ZERO;
				BigDecimal[][] Et = new BigDecimal[n][m];
				for (int i = 0; i < n; i++)
					for (int j = 0; j < m; j++)
						Et[i][j] = BigDecimal.ZERO;

				for (int sm = 0; sm < M; sm++) {
					IloCplex dual = new IloCplex();
					dual.setOut(new NullOutputStream());
					IloNumVar[] phi = dual.numVarArray(m,
							Double.NEGATIVE_INFINITY, 0.0);
					IloNumVar[] sigma = dual.numVarArray(n,
							Double.NEGATIVE_INFINITY, 0.0);
					IloNumVar[][] gamma = new IloNumVar[n][m];
					for (int i = 0; i < n; i++)
						for (int j = 0; j < m; j++)
							gamma[i][j] = dual.numVar(Double.NEGATIVE_INFINITY,
									0.0);
					IloLinearNumExpr dobj = dual.linearNumExpr();
					for (int j = 0; j < m; j++)
						dobj.addTerm(capacities[j], phi[j]);
					for (int i = 0; i < n; i++)
						dobj.addTerm(demands[sm][i], sigma[i]);
					for (int i = 0; i < n; i++) {
						for (int j = 0; j < m; j++) {
							dobj.addTerm(YValue[i][j] * capacities[j],
									gamma[i][j]);

							// constraint
							dual.addLe(dual.sum(phi[j], sigma[i], gamma[i][j]),
									-pi[i][j]);
						}
					}
					dual.addMaximize(dobj);

					if (dual.solve()) {
						subProblemObjValue = subProblemObjValue.add(BigDecimal
								.valueOf(dual.getObjValue()));
						for (int j = 0; j < m; j++)
							et = et.add(BigDecimal.valueOf(dual
									.getValue(phi[j]) * capacities[j]));
						for (int i = 0; i < n; i++)
							et = et.add(BigDecimal.valueOf(dual
									.getValue(sigma[i]) * demands[sm][i]));
						for (int i = 0; i < n; i++) {
							for (int j = 0; j < m; j++) {
								Et[i][j] = Et[i][j].add(BigDecimal.valueOf(dual
										.getValue(gamma[i][j])));
							}
						}
					}
					dual.end();
				}
				// end scenario enumeration

				et = et.divide(BigDecimal.valueOf(M), 5,
						BigDecimal.ROUND_HALF_EVEN);
				for (int j = 0; j < m; j++) {
					for (int i = 0; i < n; i++) {
						Et[i][j] = Et[i][j].multiply(
								BigDecimal.valueOf(-capacities[j])).divide(
								BigDecimal.valueOf(M), 5,
								BigDecimal.ROUND_HALF_EVEN);
					}
				}
				BigDecimal Etx = BigDecimal.ZERO;
				for (int i = 0; i < n; i++) {
					for (int j = 0; j < m; j++) {
						Etx = Etx.add(Et[i][j].multiply(BigDecimal
								.valueOf(YValue[i][j])));
					}
				}
				double wv = et.subtract(Etx).doubleValue();
				// System.out.println("@wv = " + wv);
				// System.out.println(wv + "," + thetaValue);

				if (iter > 1)
					gap = BigDecimal
							.valueOf(wv - thetaValue)
							.divide(BigDecimal.valueOf(-thetaValue), 5,
									BigDecimal.ROUND_HALF_EVEN).doubleValue();
				System.out.println("@" + iter + ": " + gap);

				// double Phi = BigDecimal
				// .valueOf(masterObjValue)
				// .subtract(BigDecimal.valueOf(thetaValue))
				// .add(subProblemObjValue.divide(BigDecimal.valueOf(M),
				// 5, BigDecimal.ROUND_HALF_EVEN)).doubleValue();
				// System.out.println("@check Phi=" + Phi);
				// if (Phi < UB)
				// UB = Phi;
				// double gap = UB - LB;
				// System.out.println("iter = " + iter + ", gap = " + gap);
				long timeTT = System.currentTimeMillis() - start;
				if (gap < 0.005 || timeTT > 3600000) {
					for (int i = 0; i < n; i++) {
						for (int j = 0; j < m; j++) {
							if (YValue[i][j] == 1)
								System.out.print("(" + i + "," + j + ") \t");
						}
						System.out.println();
					}
					long end = System.currentTimeMillis();
					double time = end - start;
					System.out.println("time consumed: " + time);
					System.out.println("iteration: " + iter);
					break;
				} else {
					iter++;
					IloLinearNumExpr optimalCut = master.linearNumExpr();
					for (int i = 0; i < n; i++)
						for (int j = 0; j < m; j++)
							optimalCut.addTerm(Et[i][j].doubleValue(), Y[i][j]);

					master.addGe(master.sum(optimalCut, theta),
							et.doubleValue());

				}
			}
			master.end();
		} catch (Exception e) {
			e.printStackTrace();
		}
	}
}
